library(coda)
library(rjags)
library(BayesRKD)
library(splines)
library(splines2)
library(stargazer)
library(rdrobust)

setwd("~/Hitler Youth Project")

y<-read.csv("Hitler Youth Master Data.csv")

attach(y)

used_data<-data.frame(numberofpunishments,monthshy,first_draft1,first_draft2,combattant,mths_east,catholic,urban,hum_cap,class,punished)

stargazer(used_data)


#
f <- numberofpunishments ~ 1
f1 <- numberofpunishments ~ first_draft1
f2 <- numberofpunishments ~ first_draft2
f3 <- numberofpunishments ~ first_draft2+first_draft2^2


##Models from the Paper###

f1<-numberofpunishments~1
f2<-numberofpunishments~first_draft2
f3<-numberofpunishments~numberofwounds

out1<-BayesRKD(formula=f1,data=y,running.variable=y$mths,fuzzy=y$monthshy,kink.point=184,degree=5,nItr=2000)
out2<-BayesRKD(formula=f2,data=y,running.variable=y$mths,fuzzy=y$monthshy,kink.point=184,degree=5,nItr=2000)
out3<-BayesRKD(formula=f3,data=y,running.variable=y$mths,fuzzy=y$monthshy,kink.point=184,degree=5,nItr=2000)

png("Treatment Effects 1.png")
plot(out1, trace = FALSE, density = TRUE, main = "Effect of Months in Hitler Youth")
dev.off()
png("Treatment Effects 2.png")
plot(out2, trace = FALSE, density = TRUE, main = "Effect of Months in Hitler Youth (Control for Draft Date)")
dev.off()
png("Treatment Effects 3.png")
plot(out3, trace = FALSE, density = TRUE, main = "Effect of Months in Hitler Youth (Control for Wounds)")
dev.off()




###Summary RKD Stats###

int_only<-summary(out1$samples[,6])
draft<-summary(out2$samples[,7])
wound<-summary(out3$samples[,7])

summarystats<-matrix(NA,6,3)
summarystats[,1]<-int_only
summarystats[,2]<-draft
summarystats[,3]<-wound

colnames(summarystats)<-c("No Controls","Controlling for Draft Date","Controlling for Wounds")

rownames(summarystats)<-c("Min","1st Quartile","Median","Mean","3rd Quartile","Maximum")

stargazer(summarystats,summary=F)


